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ABSTRACT 

The acceleration of relativistic particles due to repeated scattering across a shock 
wave remains the most attractive model for the production of energetic cosmic rays. 
This process has been analyzed extensively during the past two decades using the "two- 
fluid" model of diffusive shock acceleration. It is well known that 1, 2, or 3 distinct 
solutions for the flow structure can be found depending on the upstream parameters. 
Interestingly, in certain cases both smooth and discontinuous transitions exist for the 
same values of the upstream parameters. However, despite the fact that such multiple 
solutions to the shock structure were known to exist, the precise nature of the critical 
conditions delineating the number and character of shock transitions has remained un- 
clear, mainly due to the inappropriate choice of parameters used in the determination of 
the upstream boundary conditions. In this paper we derive the exact critical conditions 
by reformulating the upstream boundary conditions in terms of two individual Mach 
numbers defined with respect to the cosmic-ray and gas sound speeds, respectively. 
The gas and cosmic-ray adiabatic indices are assumed to remain constant throughout 
the flow, although they may have arbitrary, independent values. Our results provide 
for the first time a complete, analytical classification of the parameter space of shock 
transitions in the two-fluid model. We use our formalism to analyze the possible shock 
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structures for various values of the cosmic-ray and gas adiabatic indices. When multiple 
solutions are possible, we propose using the associated entropy distributions as a means 
for indentifying the most stable configuration. 

Subject headings: acceleration of particles — cosmic rays — methods: analytical - 
shock waves 

1. INTRODUCTION 

It is now widely accepted that acceleration in supernova-driven shock waves plays an important 
role in the production of the observed cosmic ray spectrum up to energies of ~ 10 15 eV (Heavens 
1984a; Ko 1995a), and it is plausible that acceleration in shocks near stellar-size compact objects 
can produce most of the cosmic radiation observed at higher energies (Jones & Ellison 1991). In 
the generic shock acceleration model, cosmic rays scatter elastically with magnetic irregularities 
(MHD waves) that are frozen into the background (thermal) gas (Gleeson & Axford 1967; Skilling 
1975). In crossing the shock, these waves experience the same compression and deceleration as the 
background gas, if the speed of the waves with respect to the gas (roughly the Alfven speed) is 
negligible compared with the flow velocity (Achterberg 1987). The convergence of the scattering 
centers in the shock creates a situation where the cosmic rays gain energy systematically each 
time they cross the shock. Since the cosmic rays are able to diffuse spatially, they can cross the 
shock many times. In this process, an exponentially small fraction of the cosmic rays experience an 
exponentially large increase in their momentum due to repeated shock crossings. The characteristic 
spectrum resulting from this first-order Fermi process is a power-law in momentum (Krymskii 1977; 
Bell 1978a, b; Blandford & Ostriker 1978). 

It was recognized early on that if the cosmic rays in the downstream region carry away a signif- 
icant fraction of the momentum flux supplied by the incident (upstream) gas, then the dynamical 
effect of the cosmic-ray pressure must be included in order to obtain an accurate description of the 
shock structure (Axford, Leer, & Skadron 1977). In this scenario, the coupled nonlinear problem 
of the gas dynamics and the energization of the cosmic rays must be treated in a self-consistent 
manner. A great deal of attention has been focused on the "two-fluid" model for diffusive shock 
acceleration as a possible description for the self-consistent cosmic-ray modified shock problem. In 
this steady-state theory, first analyzed in detail by Drury & Volk (1981, hereafter DV), the cosmic 
rays and the background gas are modeled as interacting fluids with constant specific heat ratios 
7 C and 7 9 , respectively. The coupling between the cosmic rays and the gas is provided by MHD 
waves, which serve as scattering centers but are otherwise ignored. The cosmic rays are treated 
as massless particles, and second-order Fermi acceleration due to the stochastic wave propagation 
is ignored. Within the context of the two-fluid model, DV were able to demonstrate the existence 
of multiple (up to 3) distinct dynamical solutions for certain upstream boundary conditions. The 
solutions include flows that are smooth everywhere as well as flows that contain discontinuous, 
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gas-mediated "subshocks." Subsequently, multiple solutions have also been obtained in modified 
two-fluid models that incorporate a source term representing the injection of seed cosmic rays (Ko, 
Chan, & Webb 1997; Zank, Webb, & Donohue 1993). Only one solution can be realized in a given 
flow, but without incorporating additional physics one cannot determine which solution it will 
be. The two-fluid model has been extended to incorporate a quantitative treatment of the MHD 
wave field by McKenzie & Volk (1982) and Volk, Drury, k McKenzie (1984) using a "three-fluid" 
approach. 

During the intervening decades, a great deal of effort has been expended on analyzing the 
structure and the stability of cosmic-ray modified shocks (see Jones & Ellison 1991 and Ko 1995b 
for reviews) . Much of this work has focused on the time-dependent behavior of the two- fluid model, 
which is known to be unstable to the development of acoustic waves (Drury & Falle 1986; Kang, 
Jones, & Ryu 1992) and magnetosonic waves (Zank, Axford, & McKenzie 1990). Ryu, Kang, 
& Jones (1993) extended the analysis of acoustic modes to include a secondary, Rayleigh- Taylor 
instability. In most cases, it is found that the cosmic-ray pressure distribution is not substantially 
modified by the instabilities. 

The two-fluid theory of DV suffers from a "closure problem" in the sense that there is not 
enough information to compute the adiabatic indices 7 S and 7 C self-consistently, and therefore they 
must be treated as free parameters (Achtcrbcrg, Blandford, & Periwal 1984; Duffy, Drury, & Volk 
1994). This has motivated the subsequent development of more complex theories that utilize a 
diffusion-convection transport equation to solve for the cosmic-ray momentum distribution along 
with the flow structure self-consistently. In these models, "seed" cosmic rays are either advected into 
the shock region from far upstream, or injected into the gas within the shock itself. Interestingly, 
Kang & Jones (1990), Achterberg, Blandford, & Periwal (1984), and Malkov (1997a; 1997b) found 
that diffusion-convection theories can also yield multiple dynamical solutions for certain values of 
the upstream parameters, in general agreement with the two-fluid model. Frank, Jones, & Ryu 
(1995) have obtained numerical solutions to the time-dependent diffusion-convection equation for 
oblique cosmic-ray modified shocks that are in agreement with the predictions of the steady-state 
two-fluid model. These studies suggest that, despite its shortcomings, the two-fluid theory remains 
one of the most useful tools available for analyzing the acceleration of cosmic rays in shocks waves 
(Ko 1995b). 

In their approach to modeling the diffusive acceleration of cosmic rays, DV stated the upstream 
boundary conditions for the incident flow in terms of the total Mach number, 

M = — , a^/^ + ^, (1.1) 
a y p p 

and the ratio of the cosmic-ray pressure to the total pressure, 

Q = y, P=P c + Pg, (1.2) 

where v, a, p, P, P c , and P g denote the flow velocity of the background gas, the total sound 
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speed, the gas density, and the total, cosmic-ray, and gas pressures, respectively. We will use the 
subscripts "0" and "1" to denote quantities associated with the far upstream and downstream 
regions, respectively. DV described the incident flow conditions by selecting values for M$ and 
Qo. Once these parameters have been specified, the determination of the flow structure (and in 
particular the number of possible solutions) in the simplest form of the two-fluid model requires 
several stages of root finding. Since the analysis is inherently numerical in nature, the results arc 
usually stated only for specific upstream conditions. 

The characterization of the upstream conditions in terms of Mq and Qo employed by DV turns 
out to be an inconvenient choice from the point of view of finding exact critical relations describing 
the number of possible flow solutions for given upstream conditions. As an alternative approach, it 
is possible to work in terms of the individual gas and cosmic-ray Mach numbers, defined respectively 
by 

M g = — , M c = - , (1.3) 



a g a c 



where 



denote the gas and cosmic-ray sound speeds, respectively. According to equation (1.1), a 2 = a 2 +a 2 , 
and therefore the Mach numbers M g and M c are related to the total Mach number M and the 
pressure ratio Q via 



M- 2 = M- 2 + M- 2 , Q = ,. Jf, ! - 9 , (1.5) 



lc M 2 + lg M 2 g' 



or, equivalently, 



m °={ i+ i^qT m - M -( i+ ^f v (i - 6) 

Since these equations apply everywhere in the flow, the boundary conditions in the two-fluid 
model can evidently be expressed by selecting values for any two of the four upstream param- 
eters (Mo, Qo, M g o, Mgo). In their work, DV described the upstream conditions using (Mo,Qo)j 
whereas Ko, Chan, & Webb (1997) and Axford, Leer, & McKenzie (1982) used (M g0 ,Qo)- Another 
alternative, which apparently has not been considered before, is to use the parameters (M g o, M c q). 
Although these choices are all equivalent from a physical point of view, we demonstrate below that 
the set (Mgo,M c o) is the most advantageous mathematically because it allows us to derive exact 
constraint curves that clearly delineate the regions of various possible behavior in the parameter 
space of the two-fluid model. This approach exploits the formal symmetry between the cosmic-ray 
quantities and the gas quantities as they appear in the expressions describing the asymptotic states 
of the flow. 

The remainder of the paper is organized as follows. In § 2 we discuss the transport equation 
for the cosmic rays and derive the associated moment equation for the variation of the cosmic-ray 
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energy density. In § 3 we employ momentum and energy conservation to obtain an exact result 
for the critical upstream cosmic-ray Mach number that determines whether smooth flow is possible 
for a given value of the upstream gas Mach number. In § 4 we derive exact critical conditions 
for the existence of multiple solutions containing a discontinuous, gas-mediated subshock. The 
resulting curves are plotted and analyzed for various values of the adiabatic indices 7 9 and 7 C . In 
§ 5 we present specific examples of multiple-solution flows that verify the predictions made using 
our analytical critical conditions. We conclude in § 6 with a general discussion of our results and 
their significance for the theory of diffusive cosmic-ray acceleration. 



2. GOVERNING EQUATIONS 

The two-fluid model is developed by treating the cosmic rays as a fluid with energy density 
comparable to that of the background gas, but possessing negligible mass. In this section we review 
the basic equations relevant for the two-fluid model. For integrity and clarity of presentation, we 
also include re-derivations of a few of the published results concerning the overall shock structure 
and the nature of the transonic flow. 



2.1. Lagrangian Equations 

The diffusive acceleration of energetic cosmic rays due to the convergence of scattering centers 
in a one-dimensional, plane-parallel flow is described by the transport equation (Skilling 1971; 1975) 

Df p df dv d ( df\ , , 



Dt 3 dp dx dx \ dx 

where p is the particle momentum, v(x, t) is the flow velocity of the background gas (taken to 
be positive in the direction of increasing x), n(x,p,t) is the spatial diffusion coefficient, and the 
operator 

D d d 

— = tt + Vtt- (2-2) 
Dt dt dx y J 

expresses the comoving (Lagrangian) time derivative in the frame of the gas. Equation (2.1) 
describes the effects of Fermi acceleration, bulk advection, and spatial diffusion on the direction- 
integrated (isotropic) cosmic-ray momentum distribution f(x,p,t), which is normalized so that the 
total number density of the cosmic rays is given by 



oo 



n c (x,t)= / 4TT P 2 fdp. (2.3) 

J 

Note that equation (2.1) neglects the second-order Fermi acceleration of the cosmic rays that occurs 
due to stochastic wave propagation, which is valid provided the Alfven speed va = BjyJ Airp is much 
less than the flow velocity v, where B is the magnetic field strength. Furthermore, equation (2.1) 
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does not include a particle collision term, and therefore it is not applicable to the background gas, 
which is assumed to have a thermal distribution. The momentum, mass, and energy conservation 
equations for the gas can be written in the comoving frame as 

Dv = _ldP Dp = _ dv R£l = ^ElRP (24) 

Dt p dx ' Dt 1 dx ' Dt 9 p Dt ' 1 ' ' 

respectively, where U g = P g /("f g ~ 1) * s the internal energy density of the gas. The expression 
for DUg/Dt in equation (2.4) implies a purely adiabatic variation of U g , and therefore it neglects 
any heating or cooling of the gas due to wave generation or damping. This adiabatic equation 
must be replaced with the appropriate Rankine-Hugoniot jump conditions at a discontinuous, 
gas-mediated subshock, should one occur in the flow. In the case of a relativistic subshock, the 
momentum conservation equation for the gas must be modified to reflect the anisotropy of the 
pressure distribution (e.g., Kirk & Webb 1988). 



2.2. Cosmic-Ray Energy Equation 

The pressure P c and the energy density U c associated with the isotropic cosmic-ray momentum 
distribution / are given by (Duffy, Drury, & Vdlk 1994) 

rOO A — pOO 

Pc(x,t)= —p 3 Vfdp, U c (x,t)= 4irp 2 efdp, (2.5) 

Jo 6 Jo 

where 



£ = ( 7 -l)mc 2 , V = ^-, T= ,/_P|L + i, (2.6) 

denote respectively the kinetic energy, the speed, and the Lorentz factor of a cosmic ray with 
momentum p and mass m. Although the lower bound of integration is formally taken to be p = 0, 
in practice the cosmic rays are highly relativistic particles, and therefore / vanishes for p < mc. If 
the distribution has the power-law form / oc p~ q , then we must have 4 < q < 5 in order to avoid 
divergence in the integrals for P c and U c (Achterberg, Blandford, & Periwal 1984), although this 
restriction can be lifted if cutoffs are imposed at high and/or low momentum (Kang & Jones 1990). 

We can obtain a conservation equation for the cosmic-ray energy density U c by operating on 
the transport equation (2.1) with J °° Air p 2 T dp, yielding 

du c TT dv d f_du c \ 

= Tc U c — I — I k —r ] , (2.7) 



Dt dx dx \ dx r 

where the mean diffusion coefficient R is defined by (Duffy, Drury, & Volk 1994) 

= f-p 2 T K (df/dx)d P 
{ ' } ~ f °°p2T(df/dx)dp ' 1 ' 
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and the cosmic-ray adiabatic index 7 C is defined by (Malkov & Volk 1996) 

4 i p 2 T f -f' 1 dp P c , 

^'» + » fw ^ ( * 

Note that in deriving equation (2.7), we have dropped an extra term that arises via integration by 
parts because it must vanish in order to obtain finite values for P c and U c . The integral expression 
in equation (2.9) indicates that 7 C must lie in the range 4/3 < 7 C < 5/3. It also demonstrates that 
7 C will evolve in response to changes in the shape of the momentum distribution /. The closure 
problem in the two-fluid model arises because / is not calculated at all, and therefore 7 C must be 
imposed rather than computed self-consistently. 



2.3. Eulerian Equations 

The conservation equations can be rewritten in standard Eulerian form as 

dp _ _dJ_ 

dt dx ' 



(2.10) 



d , . dl , 

<<"-> = -&• (2 - u > 

where the fluxes of mass, momentum, and total energy are given respectively by 

J = pv, (2.13) 

I = pv 2 + P g + P c , (2.14) 

E = X - pv 3 + v (P g + U g ) + v (P c + U c ) - R ^ . (2.15) 
The momentum and energy fluxes can be expressed in dimensionless form as 

l=^- = u + Vg + V c , (2.16) 

JVq 



c- E 1 2 . 7g r> 7c ,3 k d ( JvqVc 

«/^o 2 7 9 - 1 7c - 1 Jvq dx V 7c - 1 

where vq is the asymptotic upstream flow velocity and the dimensionless quantities u, V g , and V c 
are defined respectively by 

v P a P c , 

u=-, V g = -^-, V c ^-^. 2.18 

Vo JVo JVq 

Note that the definition of u implies that the incident flow has uq = 1. These relations can be used 
to rewrite equations (1.3) for the gas and cosmic-ray Mach numbers as 

M a 2 = ^-, Ml = ^-. (2.19) 
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2.4. The Dynamical Equation 

In this paper we shall adopt the two-fluid approximation in the form used by DV, and therefore 
we assume that the adiabatic indices j g and 7 C are each constant throughout the flow. The assump- 
tion of constant j g is probably reasonable since the background gas is expected to remain thermal 
and nonrelativistic at all locations. The assumption of constant 7 C is more problematic, since we 
expect the cosmic ray distribution to evolve throughout the flow in response to Fermi acceleration, 
but it is justifiable if the "seed" cosmic rays are already relativistic in the far upstream region. We 
also assume that a steady state prevails, so that the fluxes J, /, and E are all conserved. In this 
case the quantities V g and V c express the pressures of the two species relative to the upstream ram 
pressure of the gas po Vq = Jvq, where po is the asymptotic upstream mass density. The Eulerian 
frame in which we are working is necessarily the frame of the shock, since that is the only frame in 
which the flow can appear stationary (Becker 1998). In a steady state, the adiabatic variation of 
V g implied by equation (2.4) indicates that along any smooth section of the flow, the gas pressure 
can be calculated in terms of the velocity using 

-P 9 = V 9 ,{£) \ (2.20) 

where V gJf and denote fiducial quantities measured at an arbitrary, fixed location within the 
section of interest. According to equation (2.19), the associated variation of the gas Mach number 
along the smooth section of the flow is given by 

f u\ 1+79 

M 9 = K {-) > ( 2 - 21 ) 

where M g * denotes the gas Mach number at the fiducial location. Substituting for V g in equa- 
tion (2.16) using equation (2.20) and differentiating the result with respect to x yields the dynamical 
equation (Achterberg 1987; Ko, Chan, & Webb 1997) 

du dV c /dx . . 

— = % . 2.22 

dx M g 2 - 1 

Critical points occur where the numerator and denominator vanish simultaneously. The vanishing 
of the denominator implies that M g = 1 at the critical point, and therefore the critical point is also 
a gas sonic point (Axford, Leer, & McKenzie 1982). The vanishing of the numerator implies that 
dV c /dx = at the gas sonic point. 



2.5. Transonic Flow Structure 

We can rewrite the dimensionless momentum flux J by using equations (2.19) to substitute for 
V g and V c in equation (2.16), yielding 

I = + + (2.23) 
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Similarly, equation (2.17) for the dimensionless energy flux £ becomes 

1 R dP c 
7 C - 1 v dx 




(2.24) 



Using equation (2.23) to eliminate M c in equation (2.24) yields for the gradient of the cosmic-ray 
pressure 

= ^ ^ = G - r <) « 2 + + (r » - r ^ - s • < 2 - 25 > 

where 

Along any smooth section of the flow, g depends only on u by virtue of equation (2.21), which 
gives M g as a function of u. In the two-fluid model, the flow is assumed to become gradient-free 
asymptotically, so that dV c /dx — > in the far upstream and downstream regions (DV; Ko, Chan, 
& Webb 1997). The function g must therefore vanish as \x\ — > oo, and consequently we can express 
the values of X and £ in terms of the upstream Mach numbers M g o and M c q using 

j = l + ^L + ^o!, (2.27) 

and 

e i+iSUiQL „ 28) 

2 7<? - 1 7c - 1 

where we have also employed the boundary condition uq = 1. 

The critical nature of the dynamical equation (2.22) implies that g = at the gas sonic point. 
Hence if the flow is to pass smoothly through a gas sonic point, then g must vanish at three locations. 
It follows that one of the key questions concerning the flow structure is the determination of the 
number of points at which g = 0. We can address this question by differentiating equation (2.25) 
with respect to u, which yields 

where M = (M~ 2 + M~ 2 ) -1 / 2 is the total Mach number, and we have used the result 

dM~ 2 Mr 2 

implied by equation (2.21). For the second derivative of g we obtain 

d 2 g -l-7c-M- 2 ( 79 - 7c ) 



du 2 7 C 



(2.31) 



Since the cosmic rays have a higher average Lorentz factor than the thermal background gas, 
lg > 7c (cf. eq. 2.9), and therefore d 2 g/du 2 < 0, implying that g(u) is concave down as indicated 
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in Figure 1. Hence there are exactly two roots for u that yield g = 0, one given by the upstream 
velocity u = uq = 1 and the other given by the downstream velocity, denoted by u = u\. We 
therefore conclude that if the flow includes a gas sonic point, then the velocity at that point must 
be either no or u±. Consequently the flow cannot pass smoothly through a gas sonic point, as first 
pointed out by DV. 

The flows envisioned here are decelerating, and therefore the high-velocity root u = uq is 
associated with the incident flow. In this case, the fact that g(u) is concave-down implies that 
dg/du < in the upstream region, and therefore based on equation (2.29) we conclude that M > 1 
in the upstream region and M < 1 in the downstream region, with M = 1 at the peak of the curve 
where dg/du = 0. Hence the flow must contain a sonic transition with respect to the total sound 
speed a. In this sense, the flow is a "shock" whether or not it contains an actual discontinuity. 
For the flow to decelerate through a shock transition, the total upstream Mach number must 
therefore satisfy the condition Mo > 1. This constraint also implies that the upstream flow must 
be supersonic with respect to both the gas and cosmic-ray sound speeds (i.e., M g o > 1 and M c q > 1), 
since M ~ 2 = M~q + M~q. Furthermore, for a given value of M g o, the upstream cosmic-ray Mach 
number M c q must exceed the minimum value 



corresponding to the limit Mq = 1. The requirement that M g o > 1 forces us to conclude that if a gas 
sonic point exists in the flow, then it must be identical to the gradient-free asymptotic downstream 
state. Consequently the flow must either remain supersonic everywhere with respect to the gas, 
or it must cross a discontinuous, gas-mediated subshock. If M g > 1 everywhere, then the flow is 
completely smooth and the gas sonic point is "virtual," meaning that it exists in the parameter 
space, but it does not lie along the flow trajectory. In this case, the gas pressure evolves in a purely 
adiabatic fashion according to equation (2.20), although the total entropy of the combined system 
(gas plus cosmic rays) must increase as the flow crosses the shock, despite the fact that it is smooth. 
In § 3 we derive the critical value for the upstream cosmic-ray Mach number M c q that determines 
whether or not smooth flow is possible for a given value of the upstream gas Mach number M g Q. 



The overall structure of a cosmic-ray modified shock governed by the dynamical equation (2.22) 
can display a variety of qualitatively different behaviors, as first pointed out by DV. Depending on 
the upstream parameters, up to three distinct steady-state solutions are possible, although only one 
of these can be realized in a given situation. The possibilities include globally smooth solutions as 
well as solutions containing a discontinuous, gas-mediated subshock. Smooth flow is expected when 
the upstream cosmic-ray pressure P c q is sufficiently large since in this case cosmic ray diffusion is 
able to smooth out the discontinuity. In this section we utilize the critical nature of the dynamical 
equation to derive an analytic expression for the critical condition that determines when smooth 




(2.32) 



3. 



CRITICAL CONDITIONS FOR SMOOTH FLOW 
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flow is possible, as a function of the upstream (incident) Mach numbers M c q and M g Q. 



3.1. Critical Cosmic-Ray Mach Number 



Whether or not the flow contains a discontinuous, gas-mediated subshock, it must be smooth 
in the upstream region (preceding the subshock if one exists). We can therefore apply equation 
(2.21) for the variation of M g in the upstream region, where it is convenient to use the incident 
parameters u = 1 and M g0 as the fiducial quantities. The requirement that M g = 1 at the gas 
sonic point implies that the velocity there is given by 



u x = M 



-2/(l+7s) 



(3.1) 



which we refer to as the "critical velocity." If a sonic point exists in the flow, then u s must 
correspond to the downstream root of the equation g(u) = 0, i.e., u s = u±. Note that the value of 
u s depends only on M g o, and consequently it is independent of M c q. 

The asymptotic states of the flow are assumed to be gradient-free, and the critical conditions 
associated with the dynamical equation (2.22) imply that dV c /dx = at the gas sonic point. The 
constancy of £ and I therefore allows us to link upstream quantities to quantities at the gas sonic 
point by using equations (2.23), (2.24), (2.27), and (2.28) to write 



, 1 M cs 

lg 7c 



1 + 



M. 



-2 

go 



7 9 



+ 



M. 



-2 
cO 



7c 



and 



1 

2 + 77 



M 2 

cs 



1 



1 



M 



1 

2 + 7^ 



90 



1 



+ 



M, 



cO 



7c 



1 ' 



(3.2) 



(3.3) 



respectively, where M cs denotes the value of the cosmic-ray Mach number at the gas sonic point. 
If we substitute for u s using equation (3.1) and eliminate M cs by combining equations (3.2) and 
(3.3), we can solve for M c q to obtain an exact expression for the critical upstream cosmic-ray Mach 
number required for the existence of a sonic point in the asymptotic downstream region, 



M c0 = M cA = < 



where 



(«o-l-£K 2 + 



7c 



R o 



2 Ig-l) 



M 2 nn + 



go ^ lg-1 



(7c - 1) 



R (#o-l) M, 2 



go 



-1/2 



Ro = M, 



2/(l+7 9 ) 



(3.4) 
(3.5) 



Note that M c a is an explicit function of the upstream gas Mach number M g o . The interpretation 
is that if M c o = M c a for a given value of M g o, then the flow is everywhere supersonic with respect 
to the gas sound speed a g except in the far-downstream region, where it asymptotically approaches 
the gas sonic point. Surprisingly, this simple solution for M c a has apparently never before appeared 
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in the literature, probably due to the fact that the analytical form is lost when one works in terms 
of the alternative parameters Mo and Qo employed by DV. This can be clearly demonstrated by 
using the expressions 

M g0 = ( 1 + ^ V2 Mo , M c0 = ( 1 + 2* V2 M , (3.6) 

to substitute for M g o and M c o in equations (3.4) and (3.5) and then attempting to solve the resulting 
equation for either Qo or Mo- It is easy to convince oneself that it is not possible to express either 
of these quantities explicitly in terms of the other. In Figure 2 we depict the curve M c o = M c a in 
the (M s o,M c o) parameter space using equations (3.4) and (3.5) for various values of j g and j c . 



3.2. Smooth Flow Criterion 

We have determined that smooth flow into an asymptotic downstream gas sonic point is pos- 
sible if M c o = M c a- However, in order to obtain a complete understanding of the significance of 
the critical upstream cosmic-ray Mach number M c a, we must determine the nature of the flow 
when M c o ^ M c a- The resulting flow structure can be analyzed by perturbing around the state 
M c o = M c a by taking the derivative of the asymptotic downstream velocity u\ with respect to 
M c o, holding M g o constant. The fact that M g o is held fixed implies that the critical velocity u s 
also remains unchanged by virtue of equation (3.1). Upon differentiating, we obtain 



dM, 



Mg 



7c + 1 | (7 fl 



7c) {la 



l + uj 3 



IgUl) 



7 9 (7 9 



1)M 9 2 



1 — til 



u 



(3.7) 



which is always negative since the flow decelerates, and therefore u\ < 1. This indicates that if M c o 
is decreased from the value M c o = M c a for fixed M s o, then the downstream velocity u\ increases 
above the critical velocity u s , and therefore the flow is everywhere supersonic with respect to the 
gas sound speed. In this case there is no gas sonic point in the flow, and consequently a globally 
smooth solution is possible. Conversely, when M c o > M c a, a gas sonic point exists in the flow, and 
therefore the flow cannot be globally smooth because that would require smooth passage through 
a gas sonic point, which we have proven to be impossible. In this case, the flow must pass through 
a discontinuous, gas-mediated subshock. We conclude that globally smooth flow is possible in the 
region below each of the critical curves plotted in Figure 2. 

Note that in each case there is a critical value for M 9 o, denoted by M g A, to the right of which 
smooth flow is always possible for any value of M c o. This critical value is the solution to the 
equation 

2 2 



R A -1- 



1 



7.9 



M 2 gA + 



Ra 



la J 



7c - 



where 
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Ra ^ A#< 1+ *> , 



M 2 gA + 



R\ 



7 9 -lJ 



( 7c -l) = 0, (3.8) 



(3.9) 
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corresponding to the limit M c a —> oo in equation (3.4). We plot M g A as a function of ^ g and 7 C 
in Figure 3. When 7 S = 5/3 and 7 C = 4/3, we find that M g A = 12.28, in agreement with the 
numerical results of DV and Heavens (1984b). 



4. CRITICAL CONDITIONS FOR MULTIPLE SOLUTIONS 

DV discovered that two new discontinuous solutions become available in addition to either a 
smooth solution or another discontinuous solution when the upstream total Mach number Mq is 
sufficiently large. Subsequent authors have confirmed the existence of multiple dynamical solutions 
within the context of the two-fluid model (Achterberg, Blandford, & Periwal 1984; Axford, Leer, 
& McKenzie 1982; Kang & Jones 1990; Ko, Chan, & Webb 1997; Zank, Webb, & Donohue 1993). 
However, most of this work utilized numerical root-finding procedures and therefore it fails to 
provide much insight into the structure of the critical conditions that determine when multiple 
solutions are possible. We revisit the problem in this section by recasting the upstream boundary 
conditions using the same approach employed in § 3. In particular, we show that when the incident 
flow conditions are stated in terms of the upstream gas and cosmic-ray Mach numbers M g o and 
M c o, respectively, it is possible to obtain exact, analytical formulae for the critical curves that form 
the border of the region of multiple solutions in the parameter space. 



4.1. Post-Subshock Flow 

The existence of multiple dynamical solutions is connected with the presence in the flow of 
a discontinuous subshock mediated by the pressure of the gas. We can therefore derive critical 
criteria related to the multiple-solution phenomenon by focusing on the nature of the flow in the 
post-subshock region, assuming that a subshock exists in the flow. As we demonstrate in the 
Appendix, the energy, momentum, and particle fluxes for the cosmic rays and the background gas 
are independently conserved as the flow crosses the subshock. This implies that the quantities 
associated with the gas satisfy the usual Rankine-Hugoniot jump conditions, as pointed out by DV. 
Far downstream from the subshock, the flow must certainly relax into a gradient-free condition if a 
steady state prevails. In fact, it is possible to demonstrate that the entire post-subshock region is 
gradient-free, so that u =constant downstream from the subshock. This has already been shown by 
DV using a geometrical approach, but it can also be easily established using the following simple 
mathematical argument. First we combine the dynamical equation (2.22) with the definition of 
g(u) given by equation (2.25) to obtain the alternative form 

= ( 7c - 1) . (4.1) 

dx w k Mg - 1 y 1 



In the post-subshock gas, M g < 1 and therefore M < 1 regardless of the value of M c , since 
M~ 2 = M~ 2 +M~ 2 . It follows from equation (2.29) that dg/du > downstream from the subshock. 
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Referring to Figure 4, we wish to prove that the subshock transition must take the velocity u directly 
to the gradient-free asymptotic downstream root denoted by ui, so that g = on the immediate 
downstream side of the subshock. To develop the proof, let us suppose instead that g > in the 
post-subshock gas, corresponding to the post-subshock velocity u a > u\ in Figure 4. In this case, 
equation (4.1) implies that du/dx > in the downstream region, so that u increases along the flow 
direction, evolving away from the gradient-free root u\ in the post-subshock flow. Conversely, if 
g < in the post-subshock gas (corresponding to the velocity u\, < u\ in Figure 4), then du/dx < 
in the downstream region and consequently u decreases, again evolving away from the gradient-free 
root u\. Hence if the flow is to be stationary, then u must jump directly to the value ui, and the 
entire post-subshock region must therefore be gradient-free. This conclusion is valid within the 
context of the "standard" two-fluid model studied by DV, but Zank, Webb, & Donohue (1993) 
suggest that it may be violated in models including injection. 



4.2. Global Flow Structure 

We can derive an expression for g(u) suitable for use in the post-subshock region by employing 
equation (2.21) to eliminate M g in equation (2.25). This yields 

/ 1 \ u 1+lg u 1 ^ 13 

g(u) = r c J u 2 + r c iu + (r g - r c ) + ^ m2+ - £ > ( 4 - 2 ) 

where we have adopted the immediate post-subshock values u + and M g+ as the fiducial quantities in 
the smooth section of the flow downstream from the subshock. Since we have already established 
that the post-subshock flow is gradient-free, we can obtain an equation satisfied by the post- 
subshock velocity u + by writing 

g(u + ) = = Q - T c ) u% + T c lu + + (T g - T c ) - £ . (4.3) 

The gradient-free nature of the post-subshock flow also trivially implies 

ui = u + , M gl = M g+ , (4.4) 

where M g \ is the asymptotic downstream gas Mach number associated with the downstream velocity 
u\. Equation (4.3) can be interpreted as a relation for the immediate pre-subshock velocity n_ by 
utilizing the standard Rankine-Hugoniot jump conditions (Landau & Lifshitz 1975) 

u + _ 2 + { lg -l)M 2 g _ 2 + ( 7g -l)M g 2 _ 

u- (l 9 + 1) M*. ' M *+ 1 - lg + 2 lg M 2 g _ ' 

where the pre-subshock gas Mach number M 9 _ is related to n_ and M g o via 

M 2 g _ = M 2 g0 u 1 ^ , (4.6) 



-15- 



which follows from equation (2.21). 

By using equations (4.5) and (4.6) to eliminate u+, M g+ , and M 9 _, we can transform equa- 
tion (4.3) into a new equation for the pre-subshock velocity u_, which we write symbolically as 



Recall that the constants 2 and £ are functions of M c q and M g o by virtue of equations (2.27) and 
(2.28). In addition to satisfying the condition h = 0, acceptable roots for the pre-subshock velocity 
u_ must also exceed the critical velocity u s . This is because the flow must be supersonic with 
respect to the gas before crossing the subshock, if one exists. Equation (4.7) allows us to solve 

for the pre-subshock velocity IX clS cl function of the upstream cosmic-ray and gas Mach numbers 

M c q and M g o, respectively. In general, u- is a multi- valued function of M c q and M g o, and this 
results in the possibility of several distinct subshock solutions in certain regions of the parameter 
space. Fortunately, additional information is also available that can be utilized to calculate the 
shape of the critical curve in (M g o, M c q) space bordering the region of multiple subshock solutions. 
The nature of this information becomes clear when one examines the topology of the function 
h(u-) as depicted in Figures 5 and 6 for j g = 5/3 and j c = 4/3. We consider a sequence of 
situations with M g o held fixed and M c q gradually increasing from the minimum value M Cyia \ a given 
by equation (2.32), corresponding to the limit Mq = 1. Note that since M g o is held constant, the 
critical velocity u s also remains constant according to equation (3.1). Two qualitatively different 
behaviors are observed depending on whether or not M g o exceeds M g A, where M g A = 12.28 is the 
critical upstream gas Mach number for smooth flow calculated using equation (3.8) with j g = 5/3 
and 7 C = 4/3. 

In Figure 5a we plot ftasa function of u_ and M c q for M g o = 8, which yields for the critical 
velocity u s = 0.210. In this case the niiniinum upstream cosmic-ray Mach. number is M^c,mm — 

1.008. 

Note that initially, for small M c q, there is one root for u- corresponding to the single crossing of 
the line h = 0. Since this root does not exceed u s , a subshock solution is impossible and instead 
the flow must be globally smooth as discussed in § 3. The choice M g o = 8 satisfies the condition 
M g o < M g A, and therefore as M c q increases, the root for u_ eventually equals u s , which occurs 
when M c q = M c a = 4.25. In this case the subshock is located at the asymptotic downstream limit 
of the flow, and therefore it is identical to the gas sonic point. Equations (3.1) and (4.6) indicate 
that the pre-subshock gas Mach number M ff _ = 1 as expected. 

As we continue to increase M c q beyond the critical value M c a in Figure 5a, the root for u_ 
exceeds u s , and therefore the smooth solution is replaced by a subshock solution. We refer to this 



h(u-) = , 



(4.7) 



where 





(4.8) 
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solution as the "primary" subshock solution. Since M 9 _ > 1 in this region of the parameter space, 
the subshock plays a significant role in modifying the flow. If M c q is increased further, the primary 
root for U- increases slowly, the concave-down shape changes, and a new peak begins to emerge at 
large it_. The peak touches the line h = when M c q = 41.65, and therefore at this point a new 
subshock root appears with the value u_ = 0.944. The development of this new root can be clearly 
seen in Figure 56, where we rep lot h at a much smaller scale. As M c q continues to increase, the 
peak continues to rise, and the new «_ root bifurcates into two roots. Since the two new roots for 
u- are larger than the primary root, the corresponding pre-subshock gas Mach numbers are also 
larger, and therefore the two new subshocks are stronger than the primary subshock. We conclude 
that in this region of the (M g o, M c q) parameter space, three discontinuous subshock solutions arc 
possible, although only one can occur in a given situation. 

Another example is considered in Figure 6a, where we plot h as a function of u_ and M c q for 
M g o = 13, which yields u s = 0.146 and M c ^ m - m = 1.003. In this case, M g o > M g A, and consequently 
there is always a root for u_ below u s , indicating that globally smooth flow is possible for all 
values of M c q. Hence the "primary" subshock solution never appears in this example. However, 
for sufficiently large M c q, a peak develops in h just as in Figure 56. This peak rises with increasing 
M c q and eventually crosses the line h = at U- = 0.981 when M c q = 116, corresponding to the 
appearance of a new physically acceptable subshock root for ii- . This new root bifurcates into two 
roots as M c q is increased, which can be clearly seen in Figure 66, where h is replotted on a much 
smaller scale. It follows that in this region of (M g o,M c o) space, two discontinuous solutions are 
possible in addition to a single globally smooth solution. It is apparent from Figures 5 and 6 that 
the onset of multiple solutions is connected with the vanishing of h coupled with the additional, 
simultaneous condition 



Equations (4.7) and (4.9) can be manipulated to obtain explicit expressions for the critical up- 
stream gas and cosmic-ray Mach numbers corresponding to the onset of multiple subshock solutions 
as functions of the pre-subshock velocity it_. These critical Mach numbers are denoted by M g B and 
M c b, respectively. The logical procedure for obtaining the relations is straightforward, although 
the algebra required is somewhat tedious. The first step in the process is to solve equations (4.7) 
and (4.9) individually to derive two separate expressions for M c b- Equation (4.7) yields 




(4.9) 



which supplements equation (4.7). 



4.3. Critical Mach Numbers for Multiple Solutions 




(4.10a) 
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where 



F 1 = 4 7ff ( 7fl -l)n! 9 , 
^ = 2 79 ( 7 , - 1) [-1 - lg + ( 7s - 1) u _] u 2 y , 
F 3 = -4 7c ( lg -l)(ul 9 -1), 



F 4 = -2u 



7s 



7c n_ (ul 9 - 1) + lg ( 7g + 1) ( 



ul 9 — «_ 



+ 7c 7^ («_ - 1) (ul 3 - 2) - 7 9 7c (2 - 5 u_ + + 2 u 



7s i o„l+T9' 



^ = 7s (7p " 1) («- " 1) «- 79 [-(7 ff + 1) (7c - 1) + (1 + 7s " 3 7c + 7 S 7c) «-] • 
The solution to equation (4.9) is given by 



M cB = M gB 



( 



[G 3 + G 4 M* B + G 5 M* 



where 



G!^-27 2 nl 9 , 



9B, 



G4 = —u 



7s 



G 2 ^7 9 (7 9 -l)^_ +279 , 
G 3 = 2 797c (-2 + nl 9 ), 

~ 2 7g 7c + (7 5 + 7c - 57 9 7c + 7g + 27g 7c) u- + 7c (7s - !) n - +79 

G 5 = 7s ^^ +279 [~7c (7s -!) + (! + 7 9 - 3 7c + 7 9 7c) «-] • 



(4.10b) 
(4.10c) 
(4.10d) 



(4.10e) 
(4.10f) 

(4.11a) 

(4.11b) 
(4.11c) 
(4.11d) 

(4.11e) 
(4.11f) 



The similarity of the dependences on M g B in equations (4.10) and (4.11) for M c b suggests that 
we can derive an exact solution for M g s as a function of u_ by equating these two expressions. 
The result obtained is 



M, 



gB 



-27V- 2 1 / 3 T 5 T~ 1/3 + 2- 1 / 3 T 6 1/3 \ 1/2 
3T 2 j 



(4.12a) 



where 



71 = ( 7g - 1) n 279 j - 7 2 (1 + 7 9 ) (1 + 7c ) + 7 9 (2 + 2 lg - 7 7c + 4 797c - 7 2 7c 



+(7s + 1) [(7s " 7c) ul 3 + 7s + 7c - 5 7 9 7c + 7 9 2 + 2 7 9 2 7c] n_ | , (4.12b) 
(7c + 1) (7 9 2 - 1) - 2 ( 79 + 1) (1 + 7s " 3 7c + 7s 7c) 



T 2 ^7s(7s-1)^ +37 



+(7s - 1) (1 + 7s - 37c + 7s 7c) u 2 



(4.12d) 
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T 3 = 2 7c (1 - j 2 g ) + (5 7c - 1 - 1 9 - 5 7ff 7c + 2 7 | 7c ) «- + (7 5 - 7c - 7 9 7c + tJ) «- - (4.12c) 



Equation (4.12) can be used to evaluate the critical upstream gas Mach number M g B as a function 
of the pre-subshock velocity U-. Once M g B is determined, we can calculate the corresponding 
critical upstream cosmic-ray Mach number M c b using either equation (4.10) or equation (4.11), 
which both yield the same result. Hence equations (4.10), (4.11), and (4.12) provide a direct 
means for calculating M g B and M c b as exact functions of u_. As an example, setting j g = 5/3, 
7c = 4/3, and u_ = 0.944 yields M g B = 8 and M c b = 41.65, in agreement with the results plotted 
in Figure 5. Likewise, setting j g = 5/3, 7 C = 4/3, and ii_ = 0.981 yields M g s = 13 and M c b = 116, 
in agreement with Figure 6. 

Although our expressions for M g B and M c b simplify considerably in the special case 7 9 = 5/3, 
7 C = 4/3, we have chosen to derive results valid for general (but constant) values of 7 S and 7 C in 
order to obtain a full understanding of the sensitivity of the critical Mach numbers to variations 
in the adiabatic indices. While admittedly somewhat complex, these equations can nonetheless 
be evaluated using a hand calculator, and replace the requirement of utilizing the root-finding 
techniques employed in most previous investigations of the multiple-solution phenomenon in the 
two-fluid model. 

In Figure 7 we plot the critical curves for the occurrence of multiple solutions using various 
values of 7 9 and 7 C . This is accomplished by evaluating M g B and M c b as parametric functions 
of the pre-subshock velocity u_ using equation (4.12) along with either equation (4.10) or (4.11). 
The critical curves denote the boundary of the wedge-shaped multiple-solution region. Inside this 
region, two new subshock solutions are possible, along with either the primary subshock solution 
or the globally smooth solution. Outside this region, the flow is either described by the primary 
subshock solution or else it is globally smooth. The lower-left-hand corner of the multiple-solution 
region curves to the left and culminates in a sharp cusp. The presence of the cusp implies that 
there is a minimum value of M g o, below which multiple solutions are never possible for any value 
of M c q. Note that the multiple-solution region becomes very narrow as the values of 7 9 and 7 C 
approach each other, suggesting that the physically allowed solutions converge on a single solution 
as the thermodynamic properties of the two populations of particles (background gas and cosmic 
rays) become more similar to each other. In the limit 7 9 = 7c , multiple solutions are not allowed 
at all, and the single available solution is either smooth or discontinuous depending on the values 
of M g0 and M c0 . 

In Figure 8 we plot all of the various critical curves for the physically important case of an 
ultrarelativistic cosmic-ray distribution (7 C = 4/3) combined with a nonrelativistic background 
gas (7 g = 5/3). This is the most fully self-consistent example of the two-fluid model, since in 



T 4 = 216 7 9 7c ( 7s - 1) T| - 16 2? - 72 lg T x T 2 T 3 ul 9 , 
T 5 = -4I? -12 7g T 2 T 3 ul 3 , 
T 6 = T 4 + (T 4 2 + 4T 5 3 ) 1/2 . 



(4.12g) 



(4.12e) 
(4.12f) 
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this case we expect that the adiabatic indices will remain constant throughout the flow, as is 
assumed in the model (see eq. [2.9]). The area of overlap between the multiple-solution region 
and the smooth-solution region implies the existence of four distinctly different domains within 
the (M g o, M c q) parameter space. Within Domain I, which lies outside both the multiple-solution 
region and the smooth-solution region, the flow must be discontinuous, with exactly one (i.e., the 
primary subshock) solution available. Domain II lies inside the multiple-solution region and outside 
the smooth-solution region, and therefore within this area of the parameter space, three distinct 
subshock solutions are possible, while globally smooth flow is impossible. Domain III is formed 
by the intersection of the multiple-solution region and the smooth-solution region, and therefore 
two discontinuous subshock solutions are available as well as one globally smooth solution. Finally, 
Domain IV lies outside the multiple-solution region and within the smooth-flow region, and therefore 
one globally smooth flow solution is possible. The existence of Domain IV is consistent with our 
expectation that smooth flow will occur for sufficiently large values of the upstream cosmic-ray 
pressure P c q due to diffusion of the cosmic rays. 

If we consider a trajectory through the (M g o, M c q) parameter space that crosses into the 
multiple-solution region, then the sequence of appearance of the subshock roots for u_ depends on 
whether the boundary is crossed through the "top" or "bottom" arcs of the wedge. We illustrate 
this phenomenon for j g = 5/3 and j c = 4/3 in Figure 8, where we consider two possible paths 
approaching the point P inside the multiple-solution region, starting outside at either points Q or 
R. The two paths cross the multiple-solution boundary on different sides of the wedge. In Figure 9 
we plot h{u-) along the segment RP (with M g o = 6.5), and in Figure 10 we plot h(u-) along 
the segment QP (with M c q = 45). Note that the primary subshock root for it_ already exists at 
points Q and R since they both lie inside Domain I. When the lower section of the multiple-solution 
boundary is crossed along segment RP (M c q = M c b = 25.9 in Fig. 9), two new roots for u- appear 
at larger values than the primary root, which is the same sequence observed in Figure 5. However, 
when the upper section of the boundary is crossed along segment QP (M g o = M g B = 5.83 in 
Fig. 10), the two new roots for u_ appear at smaller values than the primary root. Hence the 
order of appearance of the subshock roots for u_ is different along each path. Despite this path 
dependence, the actual set of roots obtained at point P is the same regardless of the approach path 
taken, and therefore there is no ambiguity regarding the possible subshock solutions available at 
any point in the (M g o,M c o) parameter space. 

In Figure 11 we present summary plots that include all of the various critical curves derived 
in §§ 3 and 4 for several different values of 7 3 and j c . Note that the area of overlap between 
the multiple-solution region and the smooth-solution region observed when j g = 5/3 and 7 C = 4/3 
rapidly disappears when the difference between ^ g and -y c is reduced, due to the increasing similarity 
between the thermodynamic properties of the cosmic rays and the background gas. Ko, Chan, & 
Webb (1997) and Bulanov & Sokolov (1984) have obtained similar parameter space plots depicting 
the critical curves for smooth flow and for the onset of multiple solutions. The parameters used 
to describe the incident flow conditions are (M 9 o,Qo) m the case of Ko, Chan, & Webb (1997) 



-20- 



and (Mo,Qo) f° r Bulanov & Sokolov (1984), where Qq and Mq refer to the incident pressure 
ratio and total Mach number, respectively (see eqs. [1.1] and [1.2]). When these parameters arc 
employed directly instead of the quantities (M g o, M c q) utilized in our work, the critical curves 
must be determined by numerical root-finding. However, if desired, equations (1.5) can be used 
to transform our exact solutions for the critical curves from the (M g0 , M c0 ) space to the (M g0 , Q ) 
and (Mq,Qq) spaces. This procedure was used to obtain the results depicted in Figure 12, which 
are identical to those presented by Ko, Chan, & Webb (1997) and Bulanov & Sokolov (1984). 



5. EXAMPLES 

Our interpretation of the various critical Mach numbers derived in §§ 3 and 4 can be verified 
by performing specific calculations of the flow dynamics for a few different values of the upstream 
Mach numbers M g o and M c $. We shall focus here on cases with 7 9 = 5/3 and 7 C = 4/3, but this is 
not essential. Along smooth sections of the flow, we can calculate the velocity profile by integrating 
the dynamical equation obtained by combining equations (2.21) and (4.1), which yields 

-l 

(5.1) 
where 

/ 1 \ 7/ 1+79 7/ 1-79 

g(u) = - r c J n 2 + r c iu + (r g - r c ) - e , (5.2) 

and we have introduced the new spatial variable 

f x dx' 

T{x) = v °L W)- (5 - 3) 

The quantities M gif and are fiducial parameters measured at an arbitrary location along the 
smooth section of interest, and the constants X and £ are functions of the incident Mach numbers 
M g o and M c $ via equations (2.27) and (2.28). As discussed in § 3, if M c q < M Cj 4, then the flow is 
everywhere supersonic with respect to the gas, and therefore we can use equation (5.1) to describe 
the structure of the entire flow by setting the fiducial parameters M g * and u* equal to the asymptotic 
upstream quantities M g o and no = 1, respectively. Conversely, if M c q > M c a, then the flow must 
contain a discontinuous, gas-mediated subshock. In this case the upstream region is governed by 
equation (5.1) with M g * = M g o and u* = uq = 1, and the downstream region is governed by 
equation (5.1) with M gif = M g+ and u* = u+. The pre-subshock and post-subshock regions arc 
linked by the jump conditions given by equations (4.5), which determine M g+ and u + and yield 
g = in the entire downstream region as expressed by equation (4.3). In order to illustrate the 
various possible behaviors predicted by our critical conditions, we shall present one example from 
each of the four domains discussed in § IV (see Fig. 8) calculated using j g = 5/3 and 7 C = 4/3. 

In Figure 13 we plot the velocity profile u(t) obtained by integrating the dynamical equa- 
tion (5.1) with M g o = 4 and M c q = 4. According to Figure 8, this point lies within Domain I, 



du 



M 9* I - 
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and therefore we expect to find a single subshock solution, and no globally smooth solution. Equa- 
tion (2.25) confirms that in this case the downstream root for g(u) is u = 0.27, and therefore 
smooth flow is impossible since this does not exceed the critical velocity u s = 0.35. Analysis of 
equation (4.7) yields a single acceptable value for the pre-subshock velocity root -u_ = 0.50, with 
an associated pre-subshock gas Mach number M g _ = 1.60. Figure 13 also includes plots of the 
gas Mach number M g (r) obtained using equation (2.21), the dimensionless gas pressure V g (r) 
obtained using equation (2.20), and the dimensionless cosmic-ray pressure V c (t) obtained using 
equation (2.16). The first three quantities exhibit a jump at the subshock, whereas the cosmic-ray 
pressure is continuous since a jump in this quantity would imply an infinite energy flux. The overall 
compression ratio is 1/ui = 3.66, and the pressures increase by the factors P g i/V g o = 9.27 and 
Vci/VcQ = 9.89 between the asymptotic upstream and downstream regions. Recall that V g and V c 
express the pressures of the two species divided by the upstream ram pressure of the gas. Hence in 
this example the two species each absorb comparable fractions of the upstream ram pressure. Note 
that the increase in the cosmic-ray pressure is entirely due to the smooth part of the transition 
upstream from the discontinuous subshock, while the gas pressure experiences most of its increase 
in crossing the subshock. 

In Figure 14 we set M g o = 6 and M c q = 60, corresponding to Domain II in Figure 8. We there- 
fore expect to find three distinct, physically acceptable subshock solutions associated with these pa- 
rameter values. Analysis of equation (2.25) verifies that no smooth solutions are possible in this case 
because the downstream root u = 0.18 is less than the critical velocity u s = 0.26. Equation (4.7) 
indicates that the three acceptable roots for the pre-subshock velocity are u_ = 0.53 , 0.72 , 0.99 , 
with associated pre-subshock gas Mach numbers M g _ = 2.54 , 3.84 , 5.94 , respectively. The gas and 
cosmic-ray pressures increase by the factors V g \/V g Q = 23 , 32 , 44 , and Vd/Vco = 2126 , 1306 , 37 , 
and the overall compression ratios are l/u\ = 5.20,4.64,3.72. Note that the solution with the 
largest cosmic-ray pressure has the largest overall compression ratio and the weakest subshock. 
This is due to the fact that the cosmic-ray pressure is amplified in the smooth-flow region upstream 
from the subshock, and this region is most extended when the flow does not encounter a subshock 
until the smallest possible value of M g _. Conversely, the solution with the largest gas pressure is 
the one with the strongest subshock and the smallest overall compression ratio. 

In Figure 15 we use the upstream parameters M g o = 12.5 and M c q = 200. This point lies 
within Domain III in Figure 8, and therefore we expect to find two distinct subshock solutions 
along with one globally smooth solution. Equation (2.25) confirms that in this case a smooth 
solution is possible since the downstream root u = 0.152 exceeds the critical velocity u s = 0.150. 
Analysis of equation (4.7) yields two acceptable values for the pre-subshock velocity given by 
u_ = 0.962,0.997. The corresponding pre-subshock gas Mach numbers are M 9 _ = 11.9,12.4, re- 
spectively. In the discontinuous (subshock) solutions the gas and cosmic-ray pressures increase by 
the factors V g i/V g Q = 188, 194, and Vd/Vco = 2011 , 170, and the overall compression ratios are 
1/ui = 4.07,3.94. The subshocks in this example are strong, and therefore in the solutions con- 
taining a subshock most of the deceleration is due to the buildup of the pressure of the background 
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gas, rather than the cosmic-ray pressure. Hence both of the subshock solutions are gas-dominated. 
In the globally smooth solution, which is cosmic-ray dominated, the pressures increase by the fac- 
tors Vgi/Vgo = 23 and V c i/V c o = 40702, and the compression ratio is l/u\ = 6.57. In this case the 
cosmic rays absorb almost all of the ram pressure of the upstream gas. 

Finally, in Figure 16 we set M g o = 14 and M c q = 20. According to Figure 8, this point lies 
within Domain IV, and therefore we expect that only a single, globally-smooth solution exists for 
these upstream parameters. This prediction is verified by equation (4.7), which confirms that no 
acceptable subshock roots for u_ exists. Furthermore, equation (2.25) indicates that smooth flow 
is possible since the downstream root u = 0.15 exceeds the critical velocity u s = 0.14. Hence the 
only acceptable solution is globally smooth, with the pressure increases given by Vgi/Vgo = 23 and 
Vd/Vco = 417. The overall compression ratio is l/u\ = 6.56 and the flow is cosmic-ray dominated. 

6. DISCUSSION 

In this paper we have obtained a number of new analytical results describing the critical 
behavior of the two-fluid model for cosmic-ray modified shocks. It is well known that in this model, 
up to three distinct solutions are possible for a given set of upstream boundary conditions. The 
behaviors of the various solutions can be quite diverse, including flows that are smooth everywhere 
as well as flows that contain a discontinuous, gas-mediated subshock. The traditional approach to 
the problem of determining the types of possible solutions, employed by Ko, Chan, & Webb (1997) 
and Bulanov & Sokolov (1984), is based on stating the upstream boundary conditions in terms of 
the incident total Mach number Mq and the incident pressure ratio Qo (see eqs. [1.1] and [1.2]). In 
this approach the determination of the available solution types requires several steps of root-finding, 
and there is no possibility of obtaining analytical expressions for the critical relationships. 

The analysis presented here utilizes a fresh approach based upon a new parameterization of 
the boundary conditions in terms of the upstream gas and cosmic-ray Mach numbers M g o and M c q, 
respectively. The analytical results we have obtained in §§ 3 and 4 for the critical upstream Mach 
numbers expressed by equations (3.4), (4.11), and (4.12) provide for the first time a systematic 
classification of the entire parameter space for the two-fluid model, which remains one of the most 
powerful and practical means available for studying the problem of cosmic-ray modified shocks. 
These expressions eliminate the need for complex root-finding procedures in order to understand 
the possible flow dynamics for a given set of upstream boundary conditions, and are made possible 
by the symmetry between the gas and cosmic-ray parameters as they appear in the expressions 
describing the asymptotic upstream and downstream states of the flow. We have compared our 
quantitative results with those of Ko, Chan, & Webb (1997) and Bulanov &; Sokolov (1984), and 
they are found to be consistent. The results are valid for arbitrary (but constant) values of the 
gas and cosmic-ray adiabatic indices 7 9 and 7 C , respectively. In § 5 we have presented numerical 
examples of flow structures obtained in each of the four parameter space domains defined in Figure 8. 
These examples verify the predictions made using our new expressions for the critical Mach numbers, 
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and confirm that the largest overall compression ratios are obtained in the globally-smooth, cosmic- 
ray dominated cases. 

The existence of multiple distinct solutions for a single set of upstream boundary conditions 
demands that we include additional physics in order to determine which solution is actually realized 
in a given situation. This question has been addressed by numerous authors using various forms 
of stability analysis as well as fully time-dependent calculations. DV speculated that when three 
distinct solutions are allowed (in Domains II and III of the parameter space plotted in Fig. 8), 
the solution with the intermediate value of the cosmic-ray pressure V c will be unstable. The 
argument is based on the idea that if the cosmic-ray pressure were to increase slightly due to a 
small perturbation, then the gas would suffer additional deceleration, leading to a further increase 
in the cosmic-ray pressure. This nonlinear process would drive the flow towards the steady-state 
solution with the largest value of V c . Conversely, a small decrease in the cosmic-ray pressure would 
decrease the deceleration, leading to a smaller value for the cosmic-ray pressure. In this case the 
flow would be driven towards the steady-state solution with the smallest value of V c . Recently, 
Mond & Drury (1998) have suggested that this type of behavior may be realized as a consequence 
of a corrugational instability. Other authors (e.g., Drury &; Falle 1986; Kang, Jones, & Ryu 1992; 
Zank, Axford, & McKenzie 1990; Ryu, Kang, & Jones 1993) have argued that the globally smooth, 
cosmic-ray dominated solutions are unstable to the evolution of MHD waves in certain situations. 
Jones & Ellison (1991) suggest that even when formally stable, the smooth solutions may not be 
realizable in nature. On the other hand, Donohue, Zank, Sz Webb (1994) report time-dependent 
simulations which seem to indicate that the smooth, cosmic-ray dominated solution is indeed the 
preferred steady-state solution in certain regions of the parameter space. Hence, despite the fact 
that much effort has been expended in analyzing the stability properties of cosmic-ray modified 
shocks, there is still no clear consensus regarding which of the steady-state solutions (if any) is 
stable and therefore physically observable for an arbitrary set of upstream conditions. 

In light of the rather contradictory state of affairs regarding the stability properties of the 
various possible dynamical solutions, we propose a new form of entropy-based stability analysis. In 
this method, the entropy of the cosmic-ray distribution is calculated by first solving the transport 
equation (2.1) for the cosmic-ray distribution / and then integrating to obtain the Boltzmann 
entropy per cosmic ray, 



where k is Boltzmann's constant, n c is the cosmic-ray number density, h is Planck's constant, 
V is the system volume, and F = f/n c . According to equation (2.3), Air p 2 F(p,x) dp gives the 
probability that a randomly selected cosmic ray at location x has momentum in the range between 
p and p + dp. The cosmic-ray entropy density S c is computed using 



where the final two terms stem from the fundamental indistinguishability of the cosmic ray particles 
(of like composition), and is necessary in order to avoid the Gibbs paradox (Reif 1965). The 




+ klnV 



(6.1) 



S c = n c S c — kn c In (n c V) + kn, 



(6.2) 
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inconvenient reference to the system volume V can be removed by combining equations (6.1) and 
(6.2) to obtain the equivalent expression 



The total entropy per particle S to t for the combined gas-particle system is calculated using 



where S g and n g denote the entropy density and the number density of the background gas, respec- 
tively. One may reasonably hypothesize that the state with the largest value for S to t will be the 
preferred state in nature. This criterion may prove useful for identifying the most stable solution 
when multiple solutions are available. We plan to pursue this line of inquiry in future work. The 
results may shed new light on the structure of cosmic-ray modified shocks. 

The authors acknowledge several stimulating conversations with Frank Jones and Truong Le. 
The authors are also grateful to the anonymous referee for several useful comments. PAB also 
acknowledges support and hospitality from NASA via the NASA/ASEE Summer Faculty Fellowship 
program at Goddard Space Flight Center. 




) + kn c . 



(6.3) 



APPENDIX 



In this section we demonstrate that the cosmic-ray energy flux must be continuous across a 
velocity discontinuity (subshock), if one is present in the flow. This is turn implies that the subshock 
must be mediated entirely by the pressure of the background gas, and therefore the discontinuity 
is governed by the classical Rankine-Hugoniot jump conditions. The cosmic-ray energy flux in the 
x direction is given by 

Fc = ~ R lt + lcVUc ' (A1) 
In a steady-state, the cosmic-ray energy equation (2.7) reduces to 

dU c dv d (_ dU c \ 

dx ^ c c dx dx \ dx J ' 

which can be combined with equation (Al) to express the derivative of F c as 

dF c dU c 

-dx- = {lc - 1)v ltx-- (A3) 
Integration in the vicinity of the subshock located at x = xq yields 



,. f X0+e dF c , f x ° +e . 1N dU c J 

lim / — — dx = lim / {^ c — l)v — — dx . (A4) 
e->o ,L. n _ e dx e^o ,L. n _ e dx 



In order to avoid unphysical divergence of F c at the subshock, U c must be continuous, and therefore 
the integrand on the right-hand side of equation (A4) is no more singular than a step function. 
This implies that in the limit e — > 0, the right-hand side vanishes, leaving 

AF C = lim F c (x + e) - F c (x - e) = . (A5) 

Hence F c remains constant across the subshock. The energy, momentum, and particle fluxes of the 
gas and the cosmic rays are therefore independently conserved across the discontinuity. This allows 
us to use the standard Rankine-Hugoniot jump conditions to describe the subshock transition in 
equations (4.5). 
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FIGURE CAPTIONS 

Fig. 1. - Topology of the function g(u) given by eq. [2.25]. The quantities uq and u\ respec- 
tively denote the upstream and downstream roots for the velocity in a globally smooth solution. 

Fig. 2. - Critical upstream cosmic-ray Mach number M c a (eq. [3.4]) for smooth flow plotted 
as a function of the upstream gas Mach number M g0 in the (M g0 ,M c0 ) parameter space for (a) 
7 9 = 5/3 and various values of 7 C as indicated; (b) 7 C = 4/3 and various values of 7 9 as indicated. 
Smooth flow is not possible in the region above each curve. 

Fig. 3. - Critical upstream gas Mach number for smooth flow M g A (eq. [3.8]) plotted as a 
function of the gas and cosmic-ray adiabatic indices -f g and 7 C , respectively. When M g Q > M g A, 
smooth flow is possible for any value of M c q. 

Fig. 4. - Schematic depiction of the function g(u) (eq. [2.25]). If the flow contains a discon- 
tinuous, gas-mediated subshock, then the velocity must jump directly to the final asymptotic value 
u\ in crossing the shock. Otherwise the flow is unstable (see the discussion in the text). 

Fig. 5. - Function h(u-) (eq. [4.8]) is plotted for the parameters (a) M g o = 8, 7 9 = 5/3, 
7 C = 4/3. The value of M c q is indicated for each curve. In this example, M g o < M g A = 12.28, and 
therefore the primary subshock solution appears when the low- velocity root u_ > u s , which occurs 
when M c q > M c a = 4.25. The same function is plotted on a smaller scale in (b), where we see that 
two new subshock roots for U- appear when M c q > 41.65. 

Fig. 6. - Function h(u-) (eq. [4.8]) is plotted for the parameters (a) M g o = 13, 7 9 = 5/3, 
7 C = 4/3. The value of M c o is indicated for each curve. In this example, M g o > M g A = 12.28, and 
therefore the primary subshock solution never appears. Hence smooth flow is possible for all values 
of M c q. The same function is plotted on a smaller scale in (b), where we see that two new subshock 
roots for u_ appear when M c q > 116. 

Fig. 7. - Critical Mach numbers M g s and M c b for the onset of multiple solutions (eqs. [4.11] 
and [4.12]) are plotted as parametric functions of the pre-subshock velocity u_ in the (M g o,M c o) 
parameter space for (a) ^y g = 5/3 and the indicated values of 7 C ; (b) 7 C = 4/3 and the indicated 
values of j g . The interior of each wedge is the multiple-solution region for the associated parameters. 

Fig. 8. - Critical upstream Mach numbers for the occurrence of multiple solutions (eqs. [4.11] 
and [4.12]; solid line) and for smooth flow (eq. [3.4]; dashed line) are plotted together in the 
(M g o,M c o) parameter space for the case 7 3 = 5/3 and 7 C = 4/3. The minimum upstream cosmic- 
ray Mach number required for decelerating flow is also shown (eq. [2.32]; dotted line). There are 
four distinct domains in the parameter space as discussed in the text. 
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Fig. 9. - Function h(u-) (eq. [4.8]) is plotted for the parameters (a) M g o = 6.5, j g = 5/3, 
7 C = 4/3 along the segment RP in Fig. 8. The values of the upstream cosmic-ray Mach number 
are M c o = 45 {solid line), M c $ = 26 {dashed line), M c $ = 15 {dotted line). When M c $ > M c b = 26, 
there are three distinct subshock solutions available. In panel (6) the same function is plotted on 
a smaller scale. 

Fig. 10. - Function h{u-) (eq. [4.8]) is plotted for the parameters (a) M c q = 45, j g = 5/3, 
7 C = 4/3 along the segment QP in Fig. 8. The values of the upstream gas Mach number are 
M g0 = 6.5 {solid line), M g0 = 5.83 {dashed line), M g0 = 4.8 {dotted line). When M g0 > M gB = 5.83, 
there are three distinct subshock solutions available. In panel (6) the same function is plotted on 
a smaller scale. 

Fig. 11. - Critical upstream Mach numbers for the occurrence of multiple solutions (eqs. [4.11] 
and [4.12]; solid line) and for smooth flow (eq. [3.4]; dashed line) are plotted together in the 
{M g0 ,M c0 ) parameter space for (a) ^ g = 5/3, 7 C = 4/3; {b) j g = 5/3, 7 C = 1.35; (c) j g = 1.6, 
7c = 4/3; {d) 7 S = 1.6, 7 C = 1.35. Also indicated is the minimum value of M c q required for 
decelerating flow (eq. [2.32]; dotted line). 

Fig. 12. - Our analytical results for the critical curves generated using equations (3.4), (4.11), 
and (4.12) are combined with equations (1.5) to create corresponding curves in the alternative 
parameter spaces {Mq,Q$) and {M 9 q,Qq) employed by Bulanov & Sokolov (1984) and Ko, Chan, 
& Webb (1997), respectively. Panel (a), with j g = 5/3 and 7 C = 4/3, is identical to Fig. 4 of Bulanov 
& Sokolov (1984). Panel (6), with -y g = 2 and 7 C = 4/3, is identical to Fig. 1(a) of Ko, Chan, & 
Webb (1997). Note that these authors generated their curves using root-finding procedures. The 
interpretation of the line styles is the same as in Fig. 11. 

Fig. 13. - Numerical solutions for (a) u, (b) M g , (c) V g , and (d) V c are plotted as functions of 
r (see eq. [5.3]). The solutions were obtained by integrating the dynamical eq. [5.1] with M g o = 4 
and M c o = 4, which corresponds to Domain I in Fig. 8. In this case one discontinuous solution is 
possible, and smooth flow is impossible. 

Fig. 14. - Same as Fig. 13, except M g o = 6 and M c q = 60, which corresponds to Domain II in 
Fig. 8. In this case three distinct discontinuous solutions are possible, and smooth flow is impossible. 
The values of the pre-subshock gas Mach number are M s _ = 2.54 {solid line), M 9 _ = 3.84 {dashed 
line), M g _ = 5.94 {dotted line). 

Fig. 15. - Same as Fig. 13, except M g o = 12.5 and M c q = 200, which corresponds to Domain 
III in Fig. 8. In this case two distinct discontinuous solutions are possible in addition to one globally 
smooth solution {solid line). The values of the pre-subshock gas Mach number are M 9 _ = 11.9 
{dashed line), M g ^ = 12.4 {dotted line). 
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Fig. 16. - Same as Fig. 13, except M g o = 14 and M c $ = 20, which corresponds to Domain IV 
in Fig. 8. In this case one globally smooth solution solution is possible, and discontinuous flow is 
impossible. 
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